The heat capacity of the restricted primitive model electrolyte 
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The constant-volume heat capacity, Cv{T, p), of the restricted primitive model (RPM) electrolyte 
is considered in the vicinity of its critical point. It is demonstrated that, despite claims, recent 
simulations for finite systems do not convincingly indicate the absence of a divergence in Cv{T, p) — 
which would point to non-Ising-type criticality. The strong qualitative difference between Cv for the 
RPM and for a Lennard- Jones fluid is shown to result from the low critical density of the former. If 
one considers the theoretically preferable configurational heat-capacity density, Cv /V , the finite-size 
results for the two systems display qualitatively similar behavior on near-critical isotherms. 



The critical behavior of Coulombic systems continues 
to be subject to debate. Whereas it is generally accepted 
that the critical behavior of the gas-liquid transition 
in simple liquids belongs to the three-dimensional (3D) 
Ising universality class, the situation in ionic solutions is 
considerably more obscure. At sufficiently low tempera- 
tures, these solutions exhibit separation into two phases 
with different density, driven primarily by the Coulom- 
bic forces between the charged constituents. Experimen- 
tally, both classical (as one might guess from the long- 
range character of the ionic forces) and Ising-type crit- 
ical behavior (as might be explained by the effects of 
Debye screening) have been reported: see, e.g., Refs. 
Other possible scenarios entail a crossover from classical 
to Ising-type behavior at considerably smaller reduced 
temperatures than in simple fl^idf or even the existence 
of a different type of criticalityld'El 

In view of the significance of electrolytes and ionic sys- 
tems in many domains, a clear understanding of their 
critical behavior is of interest. It is, therefore, discon- 
certing that even for the simplest model thought to cap- 
ture the salient features of such systems, namely, the 
restricted primitive model (RPM), the universality class 
has not yet been established beyond reasonable doubt. 
The RPM consists of a mixture of hard spheres of uni- 
form diameter a, half of which carry a charge -\-q and 
half a charge —q. Its critical behavior has been analyzed 
by both analytical and numerical means. Analytically, a 
fairly satisfactory description of the critical region (ex- 
cept for the nature of the criticality) has been obtained 
from Debye-Hiickel theory supplemented by Bjerrum's 
concept of ion pairing and allowancg for the solvation of 
dipolar-ion pairs in the ionic fluid.EI However, lack of a 
sufficiently adequate formulation at the mean-field level 
has hindered the development of a renormalization-group 
treatment, see, e.g., Ref. ^. Furthermore, simulations 
have also encountered serious difficulties, not only be- 



cause of the long-range nature of the interactions, but, in 
particular, because of the low value of the critical temper- 
ature and. the resulting presence of many strongly bound 
ion pairs. B"lj The limited statistical accuracy and range in 
system sizes that have been reliably accessed have ham- 
pered detailed numerical analysis. 

This note has been Jnspired by recent work by Val- 
leau and Torrie (VT)Jj who performed numerical sim- 
ulations of the RPM using a temperature-ancUdepaity- 
scaling Monte Carlo method. Other simulationsETHllj fo- 
cused mainly on the coexistence curve below Tc- Ex- 
perimentally, observations of the coexistence curve as 
T ^ Tc~ have been revealing of universality class (with 
/3ising — 0.326 and /Jdassicai — 5) or of crossover behav- 
ior. In simulations, however, finite-size effects preclude 
the estimatioUfOf the coexistence curve close to Tc'. Wild- 
ing and BrucaHl have devised a finite-size scaling tech- 
nique for analyzing the corresponding Monte Carlo data 
which has led to fairly precise and seemingly rather re- 
liable estimates of the critical temperature, Tc, and to 
reasonable estimates of the overall ionic critical density, 
Pc- However, their technique presupposes Ising-type crit- 
icality and has not, therefore, provided any effective cri- 
teria for ruling out (or, possibly, revealing) other types 
of criticality. 

By contrast, VtA focused on the heat capacity at con- 
stant volume, Cv{T,p), in the one-phase region both as 
a function of density, p, near and on approach to crit- 
icality from above. In a classical, or van der Waals- 
type system Cy remains finite as T — > Tc+ on the 
critical isochore, p = pc, whereas in an Ising-type sys- 
tem, Cv(T,pc) diverges to infinity, albeit weakly with 
an exponent a ~ 0.109. As one passes through Tc from 
above in a classical system, Cv{T,pc) undergoes a pos- 
itive jump discontinuity and decreases smoothly there- 
after: see Fig. l|; an Ising-type fluid exhibits a \T — Tc\~°' 
singularity falling rapidly from infinity as T decreases. 
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Accordingly, VT argued that an examination of Cy {T, p) 
for the RPM for T ^ Tc and, in particular, comparison 
with simulations of a Lennard- Jones (LJ) model fluid (for 
which Ising-type or close-to-Ising-type behavior may be 
accepted) should provide an effective diagnostic of critical 
behavior. On the basis of the simulations they undertook 
and presented, VT concluded that little if any evidence 
of a rise in specific heat was present in the RPM. This 
suggested that criticality in the RPM might be classical 
in nature qe, at least, characterized by crossover rather 
close to Tc.EI 

While we acknowledge the potential value of the VT 
approach, we find, as will be explained, that we cannot 
accept the validity of their analyses or of the conclusions 
they draw. Indeed, although Ising and classical behavior 
are essentially different in the thermodynamic limit, they 
are far more difficult to distinguish in the small systems 
that are accessible to numerical simulations. Specifically, 
VT studied the heat capacity (i) along the estimated crit- 
ical isochore (for T > Tc) and (ii) along the anticipated 
critical isotherm, for a wide range of densities. In the first 
case, as mentioned, they found no signs of a divergence 
in Cv(T, pc). In the second case, no (finite-size rounded) 
peak was seen near the critical density. It is this latter 
observation that VT advance as strong evidence against 
Ising-type critical behavior in the restricted primitive 
model, since, as they illustrated, Cy (Tc, p) m a Lennard- 
Jones fluid exhibits a clear, system-size-dependent peak 
in the vicinity of p = pc- Here, we reconsider this evi- 
dence, which stands unchallenged to date, either through 
new simulations or via a reanalysis of the VT data. 

Consider, first, the heat capacity along the critical iso- 
chore. VT observe that Cv(T,pc) increases almost lin- 
early upon approach to Tc from above, with no evidence 
of a divergence. They do, however, remark that this 
might be due to the fact that all their observed tempera- 
tures lay within the regime of finite-size rounding, where 
the correlation length is restricted by the system size. We 
feel, rather, that it is the constraint t = [T — Tc)/Tc > 
that might lead to premature conclusions. In Fig. |l| we 
show the specific heat at constant volume for an infinite- 
range, van der Waals or mean-field lattice gas (in which 
all pafitieles interact equally) for a number of system 
sizes.E3^t^ The plots for density p — pc = ^Pmax rep- 
resent the behavior along the critical isochore, whereas 
the curves for p = jpc and p = jPc (the system be- 
ing symmetric around pc) illustrate the behavior along a 
noncritical isochore. As expected, the peak heights are 
lower if p 7^ but the qualitative behavior of the specific 
heat evidently persists even for relatively large deviations 
from the critical isochore. (Thus, even moderately large 
errors in the estimate of pc for the RPM should not affect 
qualitative conclusions.) 

The crucial point, however, is that (despite the absence 
of a divergence of Cy (T, pc) in the thermodynamic limit) 
the mean-field plots display pronounced size-dependent 
maxima for T < Tc- fadeed, even though these peak 
heights must saturate,t3 whereas they diverge for an 



Ising-type system with short-range interactions, the be- 
havior of small systems is qualitatively very similar in 
both cases. In particular, the specific heats of finite 
3D Ising models apd hard-core square-well fluids display 
maxima below TcM Thus, it may be difficult to distin- 
guish the two types of behavior unless one has a suffi- 
ciently large range of system sizes to allow extrapolation 
of the peak height and position. Certainly, the linear in- 
crease of Cy(r, p) for T Tc+ for a given system size, 
as VT observed for the RPM with p ~ pc^ would seem 
to convey little information regarding the nature of the 
critical behavior. This is basically a consequence of the 
fact that for finite 3D systems the specific-heat maxima 
seem invariably to occur below the true, limiting critical 
temperature. 

In order to illustrate this point more concretely, we 
have carried out high-resolution Monte Carlo simulations 
of a discretized version of the RPM.tj This model differs 
from the continuum RPM only in that the positions of the 
ions are restricted to lattice sites: the degree of discretiza- 
tion is determined by the ratio, C,, of the ion diameter a 
to the lattice spacing a. The continuuHi-|limit is recovered 
by taking Q — > oo. It has been shownOj that already for 
the small discretization parameters C = 3, 4, and 5, this 
model exhibits a liquid-vapor transition like the contin- 
uum RPM, with a coexistence curve that approaches that 
of the continuum model very closely. We have focused 
on C = 5, and carried out histogram-reweighted grand- 
canonical simulations for simple cubic lattices of sizes up 
to L = 10(7, which corresponds to {(L/a)^ = 125 000 pos- 
sible ion positions. Periodic boundary conditions were 
employed. 

The strong ion pairing at low temperatures makes 
grand-canonical simulations especially time consuming. 
However, we view canonical simulations as inherently 
dangerous owing to the important role of density fluc- 
tuations in the vicinity of the critical point. (See also 
further comments below.) Ajletailed study of these and 
related data is in progressJl3 but preliminary examina- 
tion suggests a reduced critical temperature T* ~ 0.051 
and a critical density p* ~ 0.068. (See, e.g., Refs. ||-|| for 
the standard definitions of reduced units for the RPM.) 
Our estimates agree with the suggestion of Ref. |l^ that 
the RPM with ( = 5 has a T* sfightJg iiigher than that 
of the continuum RPM {T* ~ 0.04gQl3'y), although the 
estimate of VT is T* ~ 0.052. Figure ^ shows the spe- 
cific heat of this model for p pc over a relatively wide 
temperature range around Tc- As in Fig. system-size 
dependent maxima are indeed observed at temperatures 
that clearly approach Tc from below. Their variation 
with L suggests quite strongly that Ising-type behavior 
cannot be ruled out. 

On the other hand, VT's results for the density depen- 
dence of Cy at T = Tc present, at first sight, a more 
challenging puzzle. Indeed, their data appeared to re- 
veal a remarkable difference between the RPM and a 
Lennard- Jones fluid. For the former, Cy/iVfce decreased 
monotonically with increasing density, without any ev- 
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ident marked dependence on system size. The corre- 
sponding curves VT present for a Lennard- Jones fluid, 
on the contrary, exhibit a pronounced peak in the vicin- 
ity of the critical density pc that, furthermore, increases 
with system size. 

Why is a peak for the RPM apparently absent? With- 
out doubt, any maximum or divergence related to criti- 
cality in the RPM will be rounded and shifted in finite 
systems, irrespective of the actual universality class. Fur- 
thermore, one must be prepared for strong finite-size ef- 
fects that are likely to be distorted relative to LJ-type 
model fluids in light of the shape of the iRPM coexis- 
tence curve, which is highly asymmetric.cTEl This expec- 
tation is indeed conflrmed by Fig. ^, where the continu- 
ous plots show our simulation data for the speciflc heat 
of the discretized RPM for six different system sizes at 
T* = 0.051 ~ T*. For each system size there is a clear 
maximum, at a density that appears to approach the 
critical density from below. Unfortunately, however, the 
data of VT (shown for their system size N = 192) did not 
extend to sufficiently low densities to cover these peaks. 
Also noteworthy is that their data actually agree rather 
well with ours for the discretized RPM, at a system size 
between L = 5a and L = 6cr, even though this is consid- 
erably smaller than the dimension = (N/p)^^'^ 13a 
suggested by the particle number N = 192. This dif- 
ference might be due to the C t-jS discretization; but if 
Tc really is lower for the RPM,Ej VT's data were taken 
at an isotherm somewhat above Tc- On the other hand, 
one cannot exclude the possibility that the difference re- 
flects a relative limitation of the canonical simulations 
performed by VT. The flxed particle number may, in ef- 
fect, suppress the characteristic density fluctuations, with 
a consequential relative suppression and enhanced round- 
ing of various maxima in flnite systems. This might also 
explain why VT observe negligible changes in Cy when 
increasing the number of ions by 50% from 128 to 192. 
However, the precise nature of the flnite-size effects them- 
selves is currently unclear in the framework of the T-and- 
p-scaling Monte Carlo used by VT. 

Now consider the qualitative differences observed by 
VT between the specific heats of the RPM and the LJ 
fluid. These differences also prove to be a consequence of 
the low critical density of the RPM, compared to pc for 
simple fluids; they do not reflect any signiflcant possible 
difference in the nature of the critical behavior of the two 
systems. To see this, recall that in comparing the config- 
urational heat capacity at constant volume for different 
fluids, or for a fluid and a lattice model, the more basic 
quantity for criticality and phase separation is the heat 
capacity per unit volume (or Cv density) rather than the 
speciflc heat (or heat capacity per particle)£j\ In addition 
to the arguments presented in Ref. |l^, namely the greater 
naturalness of the grand-canonical ensemble and the ex- 
pectation that spatial fluctuations most directly charac- 
terize critical divergences, note that the fleld-theoretic 
viewpoint of critical phenomena and the renormalization- 
group approachllS bear out the conclusion that the num- 



ber of degrees of freedom per unit volume plays the most 
fundamental role. 

Thus, in Fig. ^ we have plotted the heat-capacity 
density, Cy/yfee, for the RPM as a function of p for 
T* = 0.051, with a "harmless background" proportional 
to p subtracted. Instead of a monotonically decreasing 
function, we now obtain a net "energy fluctuation" that 
displays a clear maximum as a function of p, of a height 
that increases systematically with size. If we replot VT's 
data for the Lennard- Jones fluid in the analogous way, 
as in Fig. we find similar behavior: indeed, the sub- 
tracted heat-capacity density of the LJ fluid on the crit- 
ical isotherm exhibits clear finite-size maxima as a func- 
tion of the density, with a systematic size dependence 
that, as for the RPM, suggests a monotonic increase (as 
i or iV — > oo) peaking in the vicinity of the critical den- 
sity. 

In summary, the constant-volume heat capacity is a 
useful quantity in the study of the critical behavior of the 
restricted primitive model for ionic fluids. However, as 
we have demonstrated, care must be exercised before con- 
cluding that a maximum or a critical divergence is absent. 
Along both the critical isochore and the critical isotherm, 
numerical results for flnite systems exhibit clear max- 
ima, contrary to the suggestions of Ref. ^. Furthermore, 
inasfar as one observes an overall qualitative difference 
between the speciflc heat for the RPM and a fluid with 
Lennard- Jones interactions, the effect is primarily due 
to the large difference in critical densities. It has no 
relevance to possible differences in critical behavior — for 
which our current data allow few deflnitive statements. 
Conversely, if one considers the heat-capacity density, the 
two systems exhibit qualitatively rather similar behavior. 
Whether, ultimately, that reflects the same or distinct 
critical universality classes remains to be determined on 
the basis of fluctuations not only of the energy, as ob- 
served iii._tiie speciflc heat, but also of fluctuations in 
density.yiia 
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FIG. 1. The specific heat, Cv /Neks, for a grand-canonical 
lattice gas with identical interactions between all particle 
pairs, near the mean- field critical temperature . The solid 
curves derive (in order of increasing peak height) from finite 
systems of Nc = 10, 20, 40, 100, 200, 400, 1000 lattice sites, for 
a mean particle density p = pc = |pmax. The correspond- 
ing dashed curves pertain to densities p = |pc and p — \pc- 
The bold dotted curve represents the thermodynamic limit, 
Nc ^ oo, for p = Pc- The maxima of the finite-size curves 
saturate at approximately 1.657, rather than at the thermo- 
dynamic maximum |.t3 
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FIG. 2. The ionic specific heat, Cv /NkB, of the RPM with 
discretization parameter ( — 5, along the estimated critical 
isochore. Clear peaks rounded by finite size are evident be- 
low the estimated critical temperature (vertical dashed line), 
although the behavior above Tc conveys significantly less in- 
formation. 
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FIG. 3. The constant-volume specific heat on the estimated 
critical isotherm of the discretized RPM (for = 5) , compared 
with the corresponding VT data (Ref. 0). 
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FIG. 4. As in Fig. ^, but now the heat-capacity density is 
plotted and a linear "background term" Cy /Vks = 3.5p* has 
been subtracted. 
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FIG. 5. The heat-capacity density for a Lennard- Jones 
fluid, as derived from the data of Valleau and Torrie (Ref. |^ . 
As in Fig. ^, a "background" Cxr/Vk-B = 1.6/5* has been sub- 
tracted from Cy/V/cB. 
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